#ifndef SCORE_HELPERS_CL
#define SCORE_HELPERS_CL
#include <tables.cl>
#include <constants.cl>
#include <types.cl>

#pragma OPENCL EXTENSION cl_amd_printf : enable

#if 0
#define __BARRIER(type) \
  printf("\t__BARRIER hit by thread %d\n", get_local_id(0));\
  barrier(type);
#else
#define __BARRIER(type) \
  barrier(type);
#endif

__constant sampler_t s0 = CLK_NORMALIZED_COORDS_FALSE |
                          CLK_ADDRESS_CLAMP |
                          CLK_FILTER_NEAREST;
__constant sampler_t s1 = CLK_NORMALIZED_COORDS_FALSE |
                          CLK_ADDRESS_CLAMP |
                          CLK_FILTER_NEAREST;
__constant sampler_t s2 = CLK_NORMALIZED_COORDS_FALSE |
                          CLK_ADDRESS_CLAMP |
                          CLK_FILTER_NEAREST;
__constant sampler_t s3 = CLK_NORMALIZED_COORDS_FALSE |
                          CLK_ADDRESS_CLAMP |
                          CLK_FILTER_NEAREST;
__constant sampler_t s4 = CLK_NORMALIZED_COORDS_FALSE |
                          CLK_ADDRESS_CLAMP |
                          CLK_FILTER_NEAREST;
__constant sampler_t s5 = CLK_NORMALIZED_COORDS_FALSE |
                          CLK_ADDRESS_CLAMP |
                          CLK_FILTER_NEAREST;
__constant sampler_t s6 = CLK_NORMALIZED_COORDS_FALSE |
                          CLK_ADDRESS_CLAMP |
                          CLK_FILTER_NEAREST;
__constant sampler_t s7 = CLK_NORMALIZED_COORDS_FALSE |
                          CLK_ADDRESS_CLAMP |
                          CLK_FILTER_NEAREST;

float logChoose(unsigned int a, unsigned int b)
{
  return gLogFact[a] - (gLogFact[b] + gLogFact[a - b]);
}

float getScore(__local unsigned int* classCount, 
                __local MSScanInfo* scanInfo,
                __local unsigned int* peaksInScanRange,
                unsigned int peptideLength)
{
  float probability = 0.0f;
  unsigned int totalMatchedPeaks = 0;

  //Based on class counts, score the candidate
  //Compute pi(ti choose mi)
  for(unsigned int curClass = 0; 
      curClass < kNumClasses; 
      curClass++)
  {
    //Current class has 2^i * numPeaks/(2^numClasses - 1) peaks
    unsigned int curClassSize =
      (1 << curClass) * (scanInfo->mNumPeaks / ((1 << kNumClasses) - 1));

    probability += logChoose(curClassSize,
                             min(classCount[curClass], curClassSize));
  
    totalMatchedPeaks += classCount[curClass];
  }

  //Add unmatched peaks
  probability += logChoose(scanInfo->mNumHoles - scanInfo->mNumPeaks,
                           *peaksInScanRange - totalMatchedPeaks);

  //Divide by (T choose M)
  probability -= logChoose(scanInfo->mNumHoles, 
                           *peaksInScanRange);
  return -probability;
}

void insertScore(__local float* bestScoresVal,
                 __local unsigned int* bestScoresIdx,
                 float score,
                 unsigned int idx)
{
  //#pragma unroll kBestScoreListSize
  for(unsigned int curScore = 0;
      curScore < kBestScoreListSize; 
      curScore++)
  {
    //If the score is better, move shift everything less than our 
    //score down and insert our score
    if(score > bestScoresVal[curScore])
    {
      for(unsigned int i = kBestScoreListSize - 1; i > curScore; i--)
      {
        bestScoresVal[i] = bestScoresVal[i - 1];
        bestScoresIdx[i] = bestScoresIdx[i - 1];
      }

      bestScoresVal[curScore] = score;
      bestScoresIdx[curScore] = idx;

      break;
    }
  }
}

void combineScores(__local float* bestScoresLocal,
                   __local unsigned int* bestScoresIdxLocal,
                   __global float* bestScores,
                   __global unsigned int* bestScoresIdx)
{
#if kNumThreadsPerBlock > 1
  if(get_local_id(0) == 0)
  {
    for(unsigned int curScore = 0; curScore < kBestScoreListSize; curScore++)
    {
      float curMax = 0;
      unsigned int curMaxIndex = 0;
      unsigned int curMaxIdx = 0;

      //Find the max
      for(unsigned int i = 0; 
          i < kBestScoreListSize * kNumThreadsPerBlock;
          i++)
      {
        if(bestScoresLocal[i] > curMax)
        {
          curMax = bestScoresLocal[i];
          curMaxIdx = bestScoresIdxLocal[i];
          curMaxIndex = i;
        }
      }

      //Zero the max to find so we find the next largest number next time
      bestScoresLocal[curMaxIndex] = 0;

      //Write the max to the final list
      bestScores[curScore] = curMax;
      bestScoresIdx[curScore] = curMaxIdx;
    }
  }
#else
  for(unsigned int i = 0; i < kBestScoreListSize; i++)
  {
    bestScores[i] = bestScoresLocal[i];
    bestScoresIdx[i] = bestScoresIdxLocal[i];
  }
#endif
}

void findCandidates(__global unsigned int* peptideMasses,
                    unsigned int numPeptides,
                    unsigned int parentMass,
                    unsigned int curCharge,
                    __local unsigned int* firstCandidateIdx,
                    __local unsigned int* lastCandidateIdx)
{
  if(get_local_id(0) == 0)
  {
    //Search for the lower bound
    unsigned int minIndex = 0;
    unsigned int maxIndex = numPeptides - 1;

    unsigned int minErr = UINT_MAX;
    unsigned int minErrorIdx = 0;

#ifdef RELATIVE_PRECURSOR_TOLERANCE
    unsigned int searchMass = parentMass - 
                              kPrecursorTolerance * parentMass / 1000000;
#else
    unsigned int searchMass = parentMass - kPrecursorTolerance * curCharge;
#endif

    //Do a binary search for the lower bound
    do
    {
      unsigned int midPoint = (minIndex + maxIndex) / 2;

      unsigned int curMass = peptideMasses[midPoint];
      
      unsigned int err = abs((int)curMass - (int)searchMass);

      if(err <= minErr)
      {
        minErr = err;
        minErrorIdx = midPoint;
      }

      maxIndex = searchMass < curMass ? midPoint - 1 : maxIndex;
      minIndex = searchMass > curMass ? midPoint + 1 : minIndex;

    }while(minIndex <= maxIndex && minErr != 0);

    //Move the index to the first element in the tolerance
    while(peptideMasses[minErrorIdx] < searchMass)
    {
      minErrorIdx++;
    }

    while(minErrorIdx > 0 && 
          peptideMasses[minErrorIdx] == peptideMasses[minErrorIdx - 1])
    {
      minErrorIdx--;
    }

    *firstCandidateIdx = minErrorIdx;

    //Search for upper bound
    minIndex = 0;
    maxIndex = numPeptides - 1;

    minErr = UINT_MAX;
    minErrorIdx = 0;

#ifdef RELATIVE_PRECURSOR_TOLERANCE
    searchMass = parentMass + kPrecursorTolerance * parentMass / 1000000;
#else
    searchMass = parentMass + kPrecursorTolerance * curCharge;
#endif

    //Do a binary search for the upper bound
    do
    {
      unsigned int midPoint = (minIndex + maxIndex) / 2;

      unsigned int curMass = peptideMasses[midPoint];
      
      unsigned int err = abs((int)curMass - (int)searchMass);

      if(err <= minErr)
      {
        minErr = err;
        minErrorIdx = midPoint;
      }

      maxIndex = searchMass < curMass ? midPoint - 1 : maxIndex;
      minIndex = searchMass > curMass ? midPoint + 1 : minIndex;

    }while(minIndex <= maxIndex && minErr != 0);

    //Move the last candidate index to the last mass in the tolerance
    while(peptideMasses[minErrorIdx] > searchMass)
    {
      minErrorIdx--;
    }

    while(minErrorIdx < numPeptides - 1 && 
          peptideMasses[minErrorIdx] == peptideMasses[minErrorIdx + 1])
    {
      minErrorIdx++;
    }

    *lastCandidateIdx = minErrorIdx;
  }

  __BARRIER(CLK_LOCAL_MEM_FENCE);
}

unsigned int computeMZError(unsigned int a, unsigned int b)
{
#ifdef RELATIVE_TOLERANCE
  return 1000 * abs((int)a - (int)b) / a;
#else
  return abs((int)a - (int)b);
#endif
}

//Assumes data is uniformly distributed and searches for peak from there
bool findPeak(
#if defined(GLOBAL_MEMORY_SPECTRUM)
              __global const unsigned int* mz,
              __global const unsigned int* peaks,
#elif defined(LOCAL_MEMORY_SPECTRUM)
              __local const unsigned int* mz,
              __local const unsigned int* peaks,
#else
              read_only image2d_t mz,
              read_only image2d_t peaks,
#endif
              unsigned int numPeaks,
              unsigned int searchMZ,
              unsigned int* peakIntensity)
{
  int minIndex = 0;
  int maxIndex = numPeaks - 1;
  unsigned int minErr = UINT_MAX;
  unsigned int minErrorIdx = 0;

  do
  {
    unsigned int midPoint = (minIndex + maxIndex) / 2;

#if defined(LOCAL_MEMORY_SPECTRUM) || defined(GLOBAL_MEMORY_SPECTRUM)
    unsigned int curMZ = mz[midPoint];
#else
    int2 coord = { midPoint, get_global_id(1) };
    uint4 curPixel = read_imageui(mz, s0, coord);
    unsigned int curMZ = curPixel.x;
#endif    
    unsigned int err = computeMZError(curMZ, searchMZ);
    
    if(err <= minErr)
    {
      minErr = err;
      minErrorIdx = midPoint;
    }

    //Do we go left, right, or break?
    maxIndex = searchMZ < curMZ ? midPoint - 1 : maxIndex;
    minIndex = searchMZ > curMZ ? midPoint + 1 : minIndex;
  
  }while(minIndex <= maxIndex && minErr != 0);
  
#if defined(LOCAL_MEMORY_SPECTRUM) || defined (GLOBAL_MEMORY_SPECTRUM)
  *peakIntensity = peaks[minErrorIdx];
#else
  int2 coord = { minErrorIdx, get_global_id(1) };
  uint4 curPixel = read_imageui(peaks, s0, coord);
  *peakIntensity = curPixel.x;
#endif

  return minErr <= kPeakTolerance ? true : false;
}

void classifyCandidate(__local uchar* candidateString,
                       __local MSScanInfo* scanInfo,
                       unsigned int curPrecursorCharge,
                       unsigned int peptideLength,
                       __local unsigned int* classCounts,
#if defined(GLOBAL_MEMORY_SPECTRUM)
                       __global unsigned int* mz,
                       __global unsigned int* peaks,
#elif defined(LOCAL_MEMORY_SPECTRUM)
                       __local unsigned int* mz,
                       __local unsigned int* peaks,
#elif defined(IMAGE_SPECTRUM)
                       read_only image2d_t mz,
                       read_only image2d_t peaks,
#else
#error "Unknown peak packing mechanism"
#endif
                       unsigned int numPixels,
                       __local unsigned int* peaksInScanRange)
{
  unsigned int peaksInRangeReg = 0;
  unsigned int peakClass = 0;

  for(BYIon myIon = kBIon; myIon < kLastIon; myIon++)
  {
    unsigned int cumMass = myIon == kBIon ? 0 : kMassH2O;

    for(unsigned int i = 0; i < peptideLength; i++)
    {
      unsigned int curIon = myIon == kBIon ? i : (peptideLength - 1) - i;
      unsigned int curAA = (unsigned int)candidateString[curIon];

      cumMass += gMonoMass[curAA];

      //Assume we can have anywhere from a +1 to a +q charge-state on the 
      //fragment. Take the most intense one.
      unsigned int endCharge = curPrecursorCharge == 1 ? 
        1 : 
        curPrecursorCharge - 1;

      for(unsigned int curCharge = 1; 
          curCharge <= endCharge; 
          curCharge++)
      {
        unsigned int curPeak = (cumMass + kMassH * curCharge)/ curCharge;

        //If peak is out of scan range, decrement number of peaks in range
        if(curPeak < scanInfo->mMinMZ || curPeak > scanInfo->mMaxMZ)
        {
        }
        else
        {
          peaksInRangeReg++;

          //Look for a peak at the cumulative ion mass + charge state
          if(findPeak(mz,
                      peaks, 
                      scanInfo->mNumPeaks, 
                      curPeak,
                      &peakClass) == true)
          {
            classCounts[peakClass]++;
          }
        }
      }
    }
  }

  //Add the number of peaks we found to be in the scan range to the total
  *peaksInScanRange += peaksInRangeReg;
}

#endif

